/**
* \file DGVRadonTransform.h
*
* \brief This file defines the class for DRTs.
* \author Shekhar S. Chandra, 2006-2009
*
* This file is part of DGV Library.
*
* DGV is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* DGV is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with DGV.  If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef DGVRADONTRANSFORM_H
#define DGVRADONTRANSFORM_H

#include <typeinfo>

#include "DGVAliases.h"
#include "DGVFourierTransform.h"

extern "C"
{
    #include "radon.h" //FRTW Library
}

/*! \class DGVRadonTransform
    \brief Class designed to contain and encapsulate Discrete Radon Transforms (DRT).
        The transforms can be operated on Blitz Array or Matrices.

    Currently supports the following DRTs:
    - Finite Radon Transform of Matus & Flusser (1993).
    @author Shekhar S. Chandra, 2006-2009
*/
template <typename type>
class DGVRadonTransform
{
public:
	/*! \name Constructors */
	//@{
    /*! \fn DGVRadonTransform::DGVRadonTransform()
	   	\brief Standard Constructor
	   */
    DGVRadonTransform();
    /*! \fn DGVRadonTransform::~DGVRadonTransform()
	   	\brief Standard Destructor
	   */
    virtual ~DGVRadonTransform();
    //@}

	/*! \name Radon Transforms */
	//@{
	/*! \fn DGVRadonTransform::DiscreteRT(Array<type,2> &field, Array<type,2> &bins)
        \brief The Finite Radon Transform of Matus & Flusser (1993). It is assumed that the field size is prime.

        bins is automatically/internally allocated to the correct size.
    */
    void DiscreteRT(Array<type,2> &field, Array<type,2> &bins);
    /*! \fn DGVRadonTransform::iDiscreteRT(Array<type,2> &bins, Array<type,2> &result, bool norm = true)
        \brief The Inverse Finite Radon Transform of Matus & Flusser (1993). It is assumed that the field size is prime.

        bins is automatically/internally allocated to the correct size.
        If norm is true, then the result is divided by prime. It can be omitted and will normalise by default.
        */
    void iDiscreteRT(Array<type,2> &bins, Array<type,2> &result, bool norm = true);

    /*! \fn DGVRadonTransform::FastRT(Array<type,2> &field, Array<type,2> &bins)
        \brief The Fast Radon Transform of Matus & Flusser (1993). It is assumed that the field size is prime.

        These fast methods are suited only for numeric types. bins is automatically/internally allocated to the correct size.
    */
    void FastRT(Array<type,2> &field, Array<type,2> &bins);
    /*! \fn DGVRadonTransform::iFastRT(Array<type,2> &bins, Array<type,2> &result, bool norm = true)
        \brief The Inverse Fast Radon Transform of Matus & Flusser (1993). It is assumed that the field size is prime.

        These fast methods are suited only for numeric types. bins is automatically/internally allocated to the correct size.
        If norm is true, then the result is divided by prime. It can be omitted and will normalise by default.
        */
    void iFastRT(Array<type,2> &bins, Array<type,2> &result, bool norm = true);

    /*! \fn DGVRadonTransform::extractSlices(Array<complex<type>,2> &field, Array<complex<type>,2> &slices)
        \brief Extracts all the slice at discrete slopes from field. Uses the Fourier Slice Theorem.
    */
    void extractSlices(Array<complex<type>,2> &field, Array<complex<type>,2> &slices);
    void extractSlices(Array<type,2> &field, Array<type,2> &slices);
    /*! \fn DGVRadonTransform::getSlice(int m, Array<complex<type>,2> &field, Array<complex<type>,1> &slice)
        \brief Extracts the slice at discrete slope m from field. Uses the Fourier Slice Theorem.
    */
    void getSlice(int m, Array<complex<type>,2> &field, Array<complex<type>,1> &slice);
    void getSlice(int m, Array<type,2> &field, Array<type,1> &slice);
    /*! \fn DGVRadonTransform::setSlice(int m, Array<complex<type>,2> &field, Array<complex<type>,1> &slice)
        \brief Assigns a slice at discrete slope m to field. Uses the Fourier Slice Theorem.
    */
    void setSlice(int m, Array<complex<type>,2> &field, Array<complex<type>,1> &slice);
    void setSlice(int m, Array<type,2> &field, Array<type,1> &slice);
    //@}

    /*! \name I/O */
	//@{
	/*! \fn DGVRadonTransform::init()
        \brief Initializes all internal variables.
        */
	inline void init()
	{   Isum = initializer;		}
	/*! \fn DGVRadonTransform::setIsum(type sum)
        \brief Sets the \f$I_{sum}\f$ calculated in RT's (i.e. the total image sum).
        */
	inline void setIsum(type sum)
	{   Isum = sum;				}
	/*! \fn DGVRadonTransform::getIsum()
        \brief Returns the \f$I_{sum}\f$ calculated in RT's (i.e. the total image sum).
        */
	inline type getIsum()
	{   return Isum;            }
	/*! \fn DGVRadonTransform::ISumFromBins(Array<type,2> &bins)
        \brief Returns the \f$I_{sum}\f$ calculated in RT's (i.e. the total image sum).
        It is calculated from the first row of the FRT Space.
        */
	void ISumFromBins(Array<type,2> &bins);
	/*! \fn DGVRadonTransform::getInitializer()
        \brief Sets the initializer value for the RT's.
        MUST BE SET FOR CORRECT RESULTS WHEN USING TYPES OTHER THAN NUMBER TYPES (eg. string).
        */
	inline void setInitializer(type value)
	{   initializer = value;            }
	/*! \fn DGVRadonTransform::getInitializer()
        \brief Returns the initializer value for the RT's.
        MUST BE SET FOR CORRECT RESULTS WHEN USING TYPES OTHER THAN NUMBER TYPES (eg. string).
        */
	inline type getInitializer()
	{   return initializer;            }
	/*! \fn DGVRadonTransform::isTypeInteger()
        \brief Returns if the type of the RTs are integers.
        */
	bool isTypeInteger();
	//@}

protected:
    type Isum; /*!< \brief The sum of all the Field values */
    type initializer; /*!< The initialial value to which the RT's are to be conducted */
};

template <typename type>
DGVRadonTransform<type>::DGVRadonTransform()
{
    initializer = blitz::zero(initializer);
}

template <typename type>
DGVRadonTransform<type>::~DGVRadonTransform()
{

}

template <typename type>
void DGVRadonTransform<type>::DiscreteRT(Array<type,2> &field, Array<type,2> &bins)
{
	const int p = field.cols(), translate_x = 1;
    int translate_y, y, k;
    bins.resizeAndPreserve(p+1,p);
    Isum = initializer;
    bins = initializer;

    //Get Projections where j denotes projection number
    for(int j = 0; j < p; j ++)
    {
        y = 0;
        for(int m = 0; m < p; m++)
		{
			y += j;
			k = (m + translate_x)%p;
			for(translate_y = 0; translate_y < p; translate_y ++)
				bins(j,translate_y) += field(k, (translate_y + y)%p);
		}
    }
    //Get the last projection of rows
    for (int j = 0; j < p; j ++)
        for(k = 0; k < p; k ++)
        {
            bins(p,j) += field(j,k);
            Isum += field(j,k);
        }
}

template <typename type>
void DGVRadonTransform<type>::iDiscreteRT(Array<type,2> &bins, Array<type,2> &result, bool norm)
{
	const int p = bins.cols(), translate_x = 1;
    int translate_y, y, k;
    result.resizeAndPreserve(p,p);
    result = initializer;

	if(Isum == 0)
		ISumFromBins(bins);

    //Get Projections where j denotes projection number
    for(int j = 0; j < p; j ++)
    {
        y = 0;
        for(int m = 0; m < p; m++)
		{
			y += p-j;
			k = (m + translate_x)%p;
			for(translate_y = 0; translate_y < p; translate_y ++)
				result(j,translate_y) += bins(k, (translate_y + y)%p);
		}
    }

	for(int j = 0; j < p; j ++)
		for(k = 0; k < p; k++)
		{
			result(j,k) += bins(p,j);
			result(j,k) = static_cast<type>(result(j,k) - Isum);
			if(norm)
                result(j,k) /= p;
		}
}

template <typename type>
void DGVRadonTransform<type>::FastRT(Array<type,2> &field, Array<type,2> &bins)
{
	const int n = field.cols();

    Isum = initializer;

    if(n%2 == 1) ///Odd #rows then assume prime-length
    {
        bins.resizeAndPreserve(n+1,n);
        bins = initializer;

        if( isTypeInteger() )
            frt( reinterpret_cast<nttw_integer *>(field.data()), reinterpret_cast<nttw_integer *>(bins.data()), n );
        else
            frt_double( reinterpret_cast<double *>(field.data()), reinterpret_cast<double *>(bins.data()), n );
    }
    else
    {
        bins.resizeAndPreserve(n+n/2,n);
        bins = initializer;

        if( isTypeInteger() )
            frt_dyadic( reinterpret_cast<nttw_integer *>(field.data()), reinterpret_cast<nttw_integer *>(bins.data()), n );
        else
            frt_dyadic_double( reinterpret_cast<double *>(field.data()), reinterpret_cast<double *>(bins.data()), n );
    }
}

template <typename type>
void DGVRadonTransform<type>::iFastRT(Array<type,2> &bins, Array<type,2> &result, bool norm)
{
    const int n = bins.cols();

	if(Isum == 0)
		ISumFromBins(bins);

    result.resizeAndPreserve(n,n);
    result = initializer;

    if(n == bins.rows()-1) ///FRT Space has n+1 rows, then assume prime length
    {
        if( isTypeInteger() )
            ifrt( reinterpret_cast<nttw_integer *>(bins.data()), reinterpret_cast<nttw_integer *>(result.data()), n, norm );
        else
            ifrt_double( reinterpret_cast<double *>(bins.data()), reinterpret_cast<double *>(result.data()), n, norm );
    }
    else
    {
        if( isTypeInteger() )
            ifrt_dyadic( reinterpret_cast<nttw_integer *>(bins.data()), reinterpret_cast<nttw_integer *>(result.data()), n, norm );
        else
            ifrt_dyadic_double( reinterpret_cast<double *>(bins.data()), reinterpret_cast<double *>(result.data()), n, norm );
    }
}

template <typename type>
void DGVRadonTransform<type>::extractSlices(Array<complex<type>,2> &field, Array<complex<type>,2> &slices)
{
    const int p = field.cols();
    DGVFourierTransform<type,1> Fourier1D;
    Array<complex<type>,1> slice(p), projection(p);

    ///Pick slices and Inverse 1D FFT them
    slices.resize(p+1,p);
    for(int m = 0; m < p+1; m ++)
    {
        getSlice(m,field,slice);
        Fourier1D.iFFT(slice,projection);
        Fourier1D.normalise(projection);
        slices(m,Range::all()) = projection;
    }
}

template <typename type>
void DGVRadonTransform<type>::extractSlices(Array<type,2> &field, Array<type,2> &slices)
{
    const int p = field.cols();
    Array<type,1> slice(p);

    ///Pick slices and Inverse 1D FFT them
    slices.resize(p+1,p);
    for(int m = 0; m < p+1; m ++)
    {
        getSlice(m,field,slice);
        slices(m,Range::all()) = slice;
    }
}

template <typename type>
void DGVRadonTransform<type>::getSlice(int m, Array<complex<type>,2> &field, Array<complex<type>,1> &slice)
{
    const int p = slice.size();
    int index = 0;

    if(m < p && m >= 0)
    {
        slice(0) = field(0,0);
        for(int k = 1; k < p; k ++)
        {
            index = ( p-(k*m)%p )%p;
            slice(k) = field(index,k);
        }
    }
    else if(m == p)
        slice = field(Range::all(),0);
}

template <typename type>
void DGVRadonTransform<type>::getSlice(int m, Array<type,2> &field, Array<type,1> &slice)
{
    const int p = slice.size();
    int index = 0;

    if(m < p && m >= 0)
    {
        slice(0) = field(0,0);
        for(int k = 1; k < p; k ++)
        {
            index = ( p-(k*m)%p )%p;
            slice(k) = field(index,k);
        }
    }
    else if(m == p)
        slice = field(Range::all(),0);
}

template <typename type>
void DGVRadonTransform<type>::setSlice(int m, Array<complex<type>,2> &field, Array<complex<type>,1> &slice)
{
    const int p = slice.size();
    int index = 0;

    if(m < p && m >= 0)
    {
        field(0,0) += slice(0);
        for(int k = 1; k < p; k ++)
        {
            index = ( p-(k*m)%p )%p;
            field(index,k) = slice(k);
        }
    }
    else if(m == p)
        field(Range::all(),0) = slice;
}

template <typename type>
void DGVRadonTransform<type>::setSlice(int m, Array<type,2> &field, Array<type,1> &slice)
{
    const int p = slice.size();
    int index = 0;

    if(m < p && m >= 0)
    {
        field(0,0) += slice(0);
        for(int k = 1; k < p; k ++)
        {
            index = ( p-(k*m)%p )%p;
            field(index,k) = slice(k);
        }
    }
    else if(m == p)
        field(Range::all(),0) = slice;
}

template <typename type>
void DGVRadonTransform<type>::ISumFromBins(Array<type,2> &bins)
{
    Isum = initializer;
    for(int k = 0; k < bins.cols(); k ++)
        Isum += bins(0,k);
}

template <typename type>
bool DGVRadonTransform<type>::isTypeInteger()
{
    if( typeid(type).name() == "int" || typeid(type).name() == "long" || typeid(type).name() == "short" || typeid(type).name() == "char" ||
        typeid(type).name() == "long int" || typeid(type).name() == "long long" || typeid(type).name() == "long short" || typeid(type).name() == "unsigned char" ||
        typeid(type).name() == "unsigned int" || typeid(type).name() == "unsigned long" || typeid(type).name() == "unsigned short" || typeid(type).name() == "unsigned long long")
        return true;
    else
        return false;
}

#endif // DGVRADONTRANSFORM_H
